A transcriptomic based deconvolution framework for assessing differentiation stages and drug responses of AML

The diagnostic spectrum for AML patients is increasingly based on genetic abnormalities due to their prognostic and predictive value. However, information on the AML blast phenotype regarding their maturational arrest has started to regain importance due to its predictive power for drug responses. Here, we deconvolute 1350 bulk RNA-seq samples from five independent AML cohorts on a single-cell healthy BM reference and demonstrate that the morphological differentiation stages (FAB) could be faithfully reconstituted using estimated cell compositions (ECCs). Moreover, we show that the ECCs reliably predict ex-vivo drug resistances as demonstrated for Venetoclax, a BCL-2 inhibitor, resistance specifically in AML with CD14+ monocyte phenotype. We validate these predictions using LUMC proteomics data by showing that BCL-2 protein abundance is split into two distinct clusters for NPM1-mutated AML at the extremes of CD14+ monocyte percentages, which could be crucial for the Venetoclax dosing patients. Our results suggest that Venetoclax resistance predictions can also be extended to AML without recurrent genetic abnormalities and possibly to MDS-related and secondary AML. Lastly, we show that CD14+ monocytic dominated Ven/Aza treated patients have significantly lower overall survival. Collectively, we propose a framework for allowing a joint mutation and maturation stage modeling that could be used as a blueprint for testing sensitivity for new agents across the various subtypes of AML.

Published in partnership with The Hormel Institute, University of Minnesota https://doi.org/10.1038/s41698-024-00596-9 A transcriptomic based deconvolution framework for assessing differentiation stages and drug responses of AML Check for updates E. Onur Karakaslar 1,2,3 , Jeppe F. Severens 1,2,3 , Elena Sánchez-López 3,4 , Peter A. van Veelen 4 , Mihaela Zlei 5,6 , Jacques J. M. van Dongen 6,7 , Annemarie M. Otte 4 , Constantijn J. M. Halkes 8 , Peter van Balen 8 , Hendrik Veelken 3,8 , Marcel J. T. Reinders 1,2,3 , Marieke Griffioen 8 & Erik B. van den Akker 1,2,3 The diagnostic spectrum for AML patients is increasingly based on genetic abnormalities due to their prognostic and predictive value.However, information on the AML blast phenotype regarding their maturational arrest has started to regain importance due to its predictive power for drug responses.Here, we deconvolute 1350 bulk RNA-seq samples from five independent AML cohorts on a singlecell healthy BM reference and demonstrate that the morphological differentiation stages (FAB) could be faithfully reconstituted using estimated cell compositions (ECCs).Moreover, we show that the ECCs reliably predict ex-vivo drug resistances as demonstrated for Venetoclax, a BCL-2 inhibitor, resistance specifically in AML with CD14+ monocyte phenotype.We validate these predictions using LUMC proteomics data by showing that BCL-2 protein abundance is split into two distinct clusters for NPM1-mutated AML at the extremes of CD14+ monocyte percentages, which could be crucial for the Venetoclax dosing patients.Our results suggest that Venetoclax resistance predictions can also be extended to AML without recurrent genetic abnormalities and possibly to MDS-related and secondary AML.Lastly, we show that CD14+ monocytic dominated Ven/Aza treated patients have significantly lower overall survival.Collectively, we propose a framework for allowing a joint mutation and maturation stage modeling that could be used as a blueprint for testing sensitivity for new agents across the various subtypes of AML.
Acute myeloid leukemia (AML) is an aggressive hematological cancer of the myeloid lineage.AML is caused by a combination of relatively few genetic alterations that are predominantly somatically acquired and cooperatively induce a maturation arrest in combination with rapid uncontrolled proliferation of immature myeloid precursor cells.The prognosis of AML is highly dependent on the presence of such recurrent genetic alterations and varies from > 90% cure rates to <10% 1 .Therefore, the current WHO classification primarily defines AML subtypes according to the presence of eleven recurrent genetic aberrations (RGA) changes and an added heterogeneous umbrella subtype composed of a highly diverse set of relatively rare RGAs [2][3][4][5] .Only AML cases that lack any detected RGA are characterized according to their maturation stage according to the French-American-British (FAB) cytomorphological/cytochemical classification 6 .
Recently, the maturation and differentiation stage of AML blasts has gained importance due to a striking association with sensitivity or resistance to new drugs [7][8][9] .AML differentiation stage can be assessed by diagnostic flow cytometry more accurately and objectively than by cytomorphological examination alone.Current standardized flow cytometry for AML diagnosis, analyzes few carefully selected differentiation markers, sufficient for accurate immunological AML classification 10 .However, gene expression profiling by bulk RNA-sequencing (RNA-seq) analyses many more markers (several thousands) and is an attractive alternative technique as it allows both calling of genetic aberrations and estimation of cell subsets, i.e. estimated cell composition (ECC) from the same sample [11][12][13] .Reported attempts to estimate ECCs by deconvoluting bulk AML samples utilizing single-cell RNA-seq as an in-silico reference, were mostly focused on detection of survival differences without performing thorough validation of the ECCs and/or used leukemic samples as reference [14][15][16] , which prevents assessing whether ECCs are tissue-or sample-specific.
Here, we perform deconvolution, a technique in which you try to estimate each cell type within the total gene expression profile of a bulk mixture, of 1350 AML transcriptomic samples via a healthy single-cell reference while validating our findings with our in-house (LUMC) flow cytometry data.Of note, we demonstrate that the ECCs recapitulate the entire FAB landscape (M0-M7).Then, using these ECCs we predict ex-vivo drug resistance data from literature and show the agreement of these results at protein level with the help of LUMC proteomics data in AML patients, for whom we also had acquired gene expression data.To conclude, we hereby propose a transcription-based single-cell guided deconvolution framework to assess the drug effectiveness to different maturational stages of AML.We also provide our framework as a CRAN R package available at https:// github.com/eonurk/seAMLess.

Deconvolution pipeline recapitulates healthy and malignant hematopoiesis
We initially created a healthy bone marrow (BM) reference atlas from single-cell transcriptomics data (Fig. 1a).For this purpose, we integrated data from 69 101 cells covering 439 genes from three publicly available datasets (two full-transcriptome and one targeted) from two studies 17,18 (Fig. 1b, Methods).T-cells and B-cells formed separate clusters in the UMAP plot, whereas the myeloid lineage cells clustered together.To differentiate early myelopoiesis, we distinguished >3349 hematopoietic stem cells (HSC) and 1 432 erythro-myeloid progenitors (EMP), 2 939 lymphoid multipotent progenitors (LMPP), as well as 3 508 granulocytes-monocytes progenitors (GMP) (Fig. 1c).Inter-individual variability, possibly due to age differences or technical differences, mostly affected T-and B-cells and did not influence overall clustering (Fig. S1A).Homogeneous distribution of the studies on the UMAP plot showed successful integration of the different datasets (Fig. S1B).
We next performed in silico experiments to validate our deconvolution set-up.We first used the healthy BM reference to create in silico mixed bulk samples with known cell compositions.Pseudobulk profiles were simulated with one abundant cell type (80% of the cells) with the remaining cells being a mix created by random selection (Methods).Then, MuSiC 13 with the default settings was used to deconvolute the simulated pseudobulk profiles into their respective cell types (Fig. S1C).To validate on an independent dataset, we also simulated pseudobulks from 40 000 healthy bone marrow cells across 8 donors of the human cell atlas (HCA) 19 (Fig. S1D) and performed deconvolution via MuSiC.All simulated pseudobulk profiles were successfully deconvoluted into their respective cell types (matching to the most abundant cell type) apart from EMP (Fig. 1d).For this profile, the annotations were shared mostly among EMP (14%) as well as HSC (25%) and Early Erythrocytes (23%).This discrepancy might be explained by the transcriptional similarity of these cell types (Fig. 1b).
Next, we analyzed our LUMC diagnostic flow cytometry data (EuroFlow 10 panels; see Methods) for 22 AML samples with matched bulk RNA-seq data.An overview of mean fluorescence intensity (MFI) values for these samples' abnormal cells after staining with antibodies for 33 markers distributed over 7 tubes (1 Orientation + 6 AML assignment tubes) is shown in Fig. S1E (Supplementary Table S1, Fig. S7).In line with our expectations, the most abundant cell types for all samples were in the myeloid lineage for both flow cytometry analyses and estimated cell compositions (ECCs) (Fig. 1e, Supplementary Table S4).To quantify whether monocytic AML can be accurately distinguished from AML with more stem cell-like phenotypes, we plotted CD14+ monocyte percentages as determined by deconvolution against MFI values of antibodies of monocytic markers (CD11b, CD64, IREM2, and CD14) on all BM cells without gating, and observed statistically significant correlations for 3 markers (CD11b, CD64, IREM2), with CD64 being most significant (R 2 = 0.43, P < 0.001) (Fig. S1F).Also, percentages of AML cells assigned to the monocytic subset by EuroFlow panels and ECCs showed statistically significant correlations (R 2 = 0.64, P < 0.001) (Fig. 1f).
Lastly, we downloaded publicly available TARGET AML and ALL (B-ALL and T-ALL) data and deconvoluted these samples (n = 719) with the healthy BM reference (Fig. 1g, Supplementary Table S1) to show that different leukemic phenotypes could be captured by deconvolution.The ECC of the matching cell type of origin was significantly higher for the different acute leukemias (Fig. S1G), most prominent for B-ALL, confirming the ECCs' ability to capture the patients' immune phenotypes at major cell type levels.Together, these benchmarking and validation analyses demonstrated that given a healthy single cell BM transcriptomic atlas, the cell type proportions can be recapitulated faithfully via deconvolution from bulk transcriptomic leukemia cases.

Deconvolution of bulk AML transcriptomics reveals the dominating immune fraction
To investigate heterogeneity in cell composition of AML, we next applied our framework to deconvolute five independent bulk transcriptomic AML studies, i.e.TCGA-LAML 20 (n = 151), BEAT-AML 21 (n = 460), TARGET-AML 22 (n = 187), LEUCEGENE 23 (n = 452) and our cohort LUMC 11 (n = 100) totaling in 1350 samples from 1267 patients (Fig. S2A).The results are shown in Fig. 2a as a heatmap where each sample was decomposed into the 22 cell types from the healthy BM reference.We also added information on patients' clinical blast counts, FAB 6 , WHO 2016 3 , and ELN 2017 5 classes for an overarching picture of AML landscape, and we also calculated the stemness score 24 , which is a gene expression signature for patient prognosis trained on engraftment capacity of AML in immunodeficient mice (Supplementary Table S2).Each sample was assigned to the most abundant cell type as determined by deconvolution (AML phenotype) and arranged according to their fractions within each phenotype.The data showed a clear dominance of one immune cell type for most of the cases, indicating that maturational arrests and lineage skewing are leukemic properties that can be readily assessed using transcriptome sequencing.The majority of AML cases was estimated to be dominated by myeloid cells (Fig. 2a, S2A, S2B, S3A, S3D).As notable exceptions, a few pediatric AML samples from the TARGET cohort showed ECC profiles dominated by T-cells or B-cells.As previously reported these cases can nevertheless be considered as a special subgroup of pediatric AML, and specifically T-cell dominated patients were recognized for poor survival 25 .Furthermore, in line with previous reports 17 , patients with acute promyelocytic leukemia (APL), classified as AML-M3 by FAB, were correctly assigned to have a cell composition dominated by granulocyte-monocyte progenitors (GMPs).Besides APL, there are also AML cases dominated by GMP and large groups of AML assigned as monocytic AML or AML with an earlier HSC or EMP phenotype, again demonstrating that bulk transcriptomics can be used to capture the stage of arrest of AML during hematopoiesis by deconvolution.

Deconvolution of bulk AML transcriptomics agrees upon stemness score
To comprehensively visualize changes in ECCs in relation to metadata, ECCs of AML samples were visualized using a UMAP plot.As expected, samples with similar ECCs clustered, as shown after annotating the samples for their most abundant cell type (Fig. 2b-eft panel).When superimposing the deconvoluted percentages of 6 major cell types on the same UMAP plot (Fig. 2b-right panels), gradual shifts in cell composition became apparent, in particular, APL samples populated the extreme extension of the GMP cluster.Moreover, a gradient towards monocytic outgrowth and conversely to a high stemness was also observed.Comparisons with the previously published stemness score showed an overall agreement with ECCs, with patients dominated by HSCs and EMPs having a statistically significantly higher stemness score than patients dominated by CD14+ Monocytes and GMPs (Fig. S3B).Late Erythrocyte also had high stemness score, albeit with large variation.By dividing AML into cases with high and low stemness scores, we observed similar compositional changes with abundance of HSCs and EMPs in AML with high stemness scores in all cohorts.Late erythrocytes, however, were not consistently enriched in all cohorts.Furthermore, we noticed that in TARGET, T cells are abundant in AML with high  stemness scores (Fig. S3C).This raises the question whether the stemness score, which has been trained on adult AML, is useful to predict prognosis of pediatric patients particularly considering that large subtypes of adult AML such as APL and AML with mutated NPM1 are rare in pediatric AML.The consistent enrichment of HSCs and EMPs in AML with high stemness scores in all cohorts, however, further supports correct determination of the cell composition of AML by deconvoluting bulk transcriptomics.

Estimated cell composition recapitulates FAB classes
Clear associations were also observed between cell type assignments by ECC and FAB classification status (left panel of Fig. 2c).Furthermore, distinct distributions of cell type-defining gene expression profiles became evident for each FAB AML type by plotting continuous values of deconvolutions rather than categorical assignments: M0 (minimally differentiated AML) cases had high levels of the HSC-defining signature, whereas these levels decreased and EMP-and GMP signatures appeared in M1 (AML without maturation) and M2 (AML with maturation).As expected, M3 (APL) cases were almost completely dominated by the GMP signature, while M4 (Acute Myelomonocytic Leukemia) and M5 (Acute monoblastic/monocytic leukemia) samples resembled CD14+ monocyte cells and GMP.Noticeably, a few samples that were dominated by the GMP signature were assigned to AML M4 or M5 by morphological assessment (top left corner).Further inquiry in LEUCEGENE, which used more detailed FAB annotations, revealed that these AML samples were enriched for M5a (Fisher exact p = 3.08 × 10-8) and M4Eo (Fisher exact p = 2.06 × 10-6) (Fig. S3E).AML M5a are dominated by monoblastic cells 6 and AML M4Eo are characterized by myelomonocytic marrow infiltration with eosinophils containing abnormal immature granules 26 .Both AML subtypes are more differentiated than GMP, but less differentiated and clearly distinct from AML M4 and M5b, which contain more mature promonocytic or monocytic cells.The small groups of M6 (Acute Erythroleukemia) and M7 (Acute Megakaryoblastic Leukemia) AML were dominated by the signatures of Late Erythrocytes and Megakaryocyte Progenitors, respectively.These data demonstrated and confirmed that our single cell guided deconvolution strategy successfully captures the maturational arrest of AML cells at different differentiation stages of hematopoiesis.

Estimated cell composition captures genetic subtype-specific resistances to various drugs
To explore whether ECCs convey information on drug resistances, ex vivo drug response data of 122 small molecule inhibitors provided as area under the curves (AUC) for 363 AML samples from BEAT were downloaded.A higher AUC indicates that cancer cells are relatively resistant since higher drug concentrations are needed to induce cell death.To understand whether drug resistance of AML samples can be predicted by ECCs, we trained random forest (RF) models per drug via a leave-one-out cross validation (LOOCV) setting.For each RF model Spearman ρ values were calculated (Supplementary Table S5-S7).The strongest correlation was observed for the BCL2 inhibitor Venetoclax 27 (ABT-199) drug (Spearman ρ = 0.509).Drugs like EGF-R inhibitor Erlotinib (Spearman ρ = 0.376) and the mTOR pathway inhibitor Rapamycin (Spearman ρ = 0.368) are amongst the top 10 drugs for which resistance could be best predicted (Fig. S4A).
Next, we asked how changes in each cell type affect drug responses and therefore univariately associated ECC to drug resistance (Fig. 3a, Supplementary Table S8).This analysis revealed that most change in resistance occurs across maturational states.For instance, CD14+ Monocytic AML are more resistant to Venetoclax, whereas more immature cell subsets are more sensitive.This trend was also clear after overlaying the AUC values of Venetoclax onto the UMAP plot (Fig. 3b) or when comparing AUC values across states (Fig. S5A) (one-way ANOVA p = 6.31 × 10 −13 ; Supplementary Table S9).

Estimated CD14+ Monocyte percentages predict Venetoclax resistance better than BCL-2 mRNA expression
Since Venetoclax is targeting the anti-apoptotic BCL-2 protein, we next checked BCL-2 mRNA expression levels in AML and overlaid the CD14+ Monocyte percentages for these cases (Fig. 3d).The data showed that low BCL2 expression indeed correlated with strong resistance to Venetoclax.However, a few samples were resistant despite relatively high BCL-2 expression (samples at upper right corner in Fig. 3d).As the majority of these cases had a CD14+ Monocyte phenotype, we univariately associated BCL-2 expression and CD14+ Monocyte percentage to investigate which factor best explains the Venetoclax response (Fig. 3e).We also associated metadata such as reported blast percentages, ELN status and primary diagnosis to determine their effect sizes and significance (Supplementary Table S12).This analysis demonstrated that both CD14+ Monocyte percentages (adjusted P = 1.7 × 10 −18 ) and BCL-2 mRNA expression (adjusted P = 4.8 × 10 −14 ) were associated significantly with the Venetoclax response.Furthermore, we investigated whether previously reported MAC Score 28 Fig. 1 | Schematic of the study and compositional validation.a The overview of the study.Integrated cells for healthy bone marrow (BM) dataset were collected from two studies 17,18 .Five independent bulk transcriptomic AML cohort (TCGA-LAML, BEAT-AML, TARGET-AML, LEUCEGENE and LUMC) samples (n = 1,350) were deconvoluted.Drug resistance data (n = 122) from BEAT-AML were predicted using ECCs.Additional ALL samples (n = 532) from TARGET study were used for validating deconvolution framework.Proteomics (n = 39) and flow cytometry samples (n = 22) from LUMC cohort were used for further validation of the framework.b UMAP plot of integrated healthy BM (n = 69,101) cells.Annotations are lifted via Azimuth framework.T and B cell subset were merged, resulting into 22 cell types (abbreviations; GMP: Granulocyte Monocyte Progenitors, LMPP: Lymphoid Primed Multipotent Progenitors, NK: Natural Killer, EMP: Erythroid Megakaryocyte Progenitor, pDC: Plasmacytoid Dendritic Cell, CLP: Common Lymphoid Progenitor, HSC: Hematopoietic Stem Cell, BaEoMA: Basophil Eosinophil Mast Progenitor, Prog Mk: Progenitor Megakaryocyte, pre-mDC: Precursor Myeloid Dendritic Cell, ASDC: AXL+ Dendritic Cell).c Barplot for cells per cell type of the healthy BM.X-axis is in log10 scale.d Heatmap showing the fraction of deconvoluted cell types for simulated pseudobulks from HCA subset.Each row indicates a simulated pseudobulk with an overabundant cell type (80%) and adds up to 1, and each column is a cell type from the healthy BM reference.e Cellular composition of the same samples via flow gating of orientation tube (ALOT) and ECCs (RNA-seq).f Percentage of monocytic subsets measured by flow cytometry (fresh material) and ECCs (thawed for RNA-seq).g Ternary plot showing the ECCs of TARGET AML and ALL cohorts (n = 719), each dot represents a sample and corners indicate a major cell type and colors indicate the primary diagnosis of each leukemic sample (Supplementary Table S1).(albeit in gene expression) was more associated to Venetoclax resistance compared to BCL2 gene expression alone (R = -0.54,and R = -0.53),but the results were similar, and the estimated CD14 Monocytes percentages had a stronger correlation (R = 0.60) (Fig. S6E).Next, we created a multivariate model to investigate whether the significance of CD14+ Monocyte percentages diminishes along with the presence of BCL-2 expression in the same model (Fig. 3f, Multivariate Venetoclax Tab in Supplementary Table 13).In this model, CD14+ Monocyte percentages remained significantly associated with Venetoclax resistance (P = 1.92 × 10 −5 ), while BCL-2 expression was below the significance threshold of 0.01 (P = 0.015).In conclusion, these results indicate that cellular composition is a more robust marker than BCL-2 mRNA expression to predict Venetoclax resistance, specifically for AML from NOS and NPM1 mutated patients.

Estimated CD14+ Monocyte percentages associates with BCL-2 protein abundance
Next, we compared the effects of BCL-2 gene expression and CD14+ Monocyte percentages with BCL-2 protein abundance within each sample.
For this purpose, we used our LUMC produced proteomics data, and after batch correction (Fig. S6A, S6B, see Methods) for matched AML cases (n = 39; LUMC) to correlate the abundance of Apoptosis regulator BCL-2 protein vs gene expression (Fig. S6C) and CD14+ Monocyte percentages (Fig. S6D).The data showed that BCL-2 gene expression and CD14+ Monocyte percentages both correlated with BCL-2 protein abundance to a similar extent (R 2 = 0.45, P < 0.001).We also overlaid plots for BCL-2 expression and CD14+ Monocyte percentages with information on patients' genetic abnormalities (Fig. 4a, b).Based on their low AUC for Venetoclax response (Fig. 3b), two AML cases with PML-RARA had high BCL-2 protein expression (n = 2, dark green).We also confirmed the strong variability in AUC for Venetoclax response within the group of AML with mutated NPM1 (n = 18, green) by showing two distinct groups at the extremes of CD14+ Monocyte percentages which correlated with BCL-2 protein expression (Fig. 4b).Within the group of AML patients with mutated NPM1, CD14+ Monocyte percentage associated stronger with both BCL-2 protein abundance (R 2 = 0.56, P < 0.001) than BCL-2 gene expression (R 2 = 0.52) (Fig. 4c, d).In conclusion, the data demonstrate the relevance of our deconvolution approach on bulk RNA-Seq to separate AML, especially those with mutated NPM1, with high and low monocyte percentages to predict the patient's response to Venetoclax.

Estimated CD14+ Monocyte phenotype captures Venetoclax treated patient prognosis
To test whether deconvolution approach associates with clinical outcome, we utilized BEAT-AML study and selected Venetoclax/Azacitidine (Ven/ Aza) treated patients with overall survival and RNA-seq data (n = 20).Deconvolution revealed that 5 of these patients had a CD14+ monocytic phenotype and these patients had a significantly worse overall survival (p = 0.047) compared to others (Fig. 4e; Supplementary Table S14).On the contrary, stemness score 24 (p = 0.11) and ELN classification (p = 0.14) were not able to fully split these patients with the notable exception of ELN favorable patients.In summary, this analysis highlights the clinical significance of employing ECC-based annotation for Ven/Aza treated patients.

Discussion
In this work, we utilized single cell guided deconvolution to decompose bulk transcriptomics data from five independent AML cohorts and show that the obtained estimated cell compositions (ECCs) faithfully reconstitute the FAB landscape (M0-M7).Moreover, using same-sample flow cytometry, we were able to validate our deconvolution framework.Hence, following our previous work using deep transcriptomics to call various types of leukemia-defining genetic aberrations 11,29 , our current findings further underpin the power of transcriptomic-based approaches as a comprehensive and versatile platform for AML diagnosis.We next illustrated the potential use of our deconvolution framework for precision medicine applications, by correlating the estimated ECC to the results of an ex vivo drug resistance screening of 122 small-molecule inhibitors in the BEAT-AML study.For the BCL-2 inhibitor Venetoclax we show that higher levels of the estimated ECC subset 'CD14+ Monocyte' correspond to a higher resistance, and intriguingly, that estimated CD14+ Monocyte levels is a better explanatory variable of resistance to Venetoclax than BCL-2 expression alone.Nevertheless, using same-sample LUMC proteomics data in 39 patients, we show that the estimated CD14+ Monocyte levels accurately mark BCL protein expression, and that for NPM1-mutated patients the presence or absence of a CD14 monocytic outgrowth corresponds with a distinct NPM1 protein abundance.Lastly, we show that CD14 monocyte phenotype correlates with poor survival outcome.Our findings may potentially have important implications on drug use especially for genetically uncharacterized patients (AML, NOS) currently accounting for ~40% of all AML as well as other well-characterized patients such as NPM1 mutated samples.Kuusanmaki et al. reported that monocytic differentiation of AML reduced sensitivity to Venetoclax ex vivo 8 , and also with a recent paper they show ex-vivo drug responses correlate with AML response in clinic 30 .Similarly, Pei et al. reported that the different monocytic subclones in vivo created resistance to Venetoclax treatment 9 , and also recently, White et al. showed that BCL-2 inhibitor resistance could be predicted via the genes associated with monocytes.Recently, using the flow cytometry data from untreated MDS patients Ganan-Gomez et al. showed that after relapse and becoming secondary AML (sAML) patients, those patients with less maturated cell types (EMP) before treatment had a faster complete remission (CR) and longer relapse free status compared to more matured cell type (GMP) with Venetoclax treatment, supporting our hypothesis that ECCs could also predict Venetoclax resistance for MDS and sAML patients 31 .Collectively, these studies suggest that Venetoclax has different resistances at different maturational stages, and especially higher resistance for patients with CD14+ Monocyte dominated phenotype as we have shown in this manuscript, and we provide an open-source framework, seAMLess, for replicating our results or applying it to other clinically relevant datasets.
A limitation of our deconvolution strategy is that it cannot distinguish the cancerous cell types as it uses a healthy bone marrow as a reference.Although it is conceptually appealing, we have three rationales behind not using a cancerous reference.First, without mutation calling for all cells, one cannot be sure whether a cell is cancerous or not.Strategies like predicting cancer cells based on their transcriptional similarity of cells with mutation calling, proposed by Van Galen et al. 17 , adds another level of ambiguity to already not perfect deconvolution pipelines.Secondly and more importantly, heterogeneity of AML causes further sub-clustering within individual AML cases (e.g.UMAP plots of Triana et al. 18 ), therefore creating a not wellcharacterized cell type signature but rather patient specific clusters 10,18 .Lastly, using healthy subsets as reference allows our framework to provide more interpretable and intuitive results for clinicians and doctors, as it reports immune phenotypes and percentages on contrary to score-based prognostic values 16,24,29 .To summarize, we believe our proposed pipeline could be a blueprint for assessing new drugs' resistances on different cell types of AML and along with our framework, they may provide better insights for clinicians and help paving the way into precision medicine in AML.

Methods
Creating the healthy BM single-cell reference We downloaded three different healthy BM datasets from two different studies, namely Van Galen et al. 17 (full-transcriptome, n = 6915), and Triana et al. 18 (full-transcriptome, n = 13,165; 462 targeted mRNA, n = 49,057).Then, all cell labels were uplifted up via Seurat package 32 (v4.0.3) default query annotation pipeline to match with the Triana's full-transcriptome cell labels as it had the most recent and detailed labels.Next, Seurat's integration pipeline with CCA was run and the cells that were labeled as doublets/ multiplets were removed from the down-stream analyses, and this yielded a healthy BM atlas of 69,130 cells in total, covering 439 genes.Also, we used the 40,000 cell subset of HCA provided by SeuratObject package (v4.0.2) for the validation analyses (Fig. 1c, d).

Different schemes of creating pseudobulks
To create a cancerous-like pseudobulk profile from healthy BM reference and HCA subset, we selected a total of 1000 cells for each profile, majority (80%) of them coming from a cell-type (over-abundant) and the rest of the cells were distributed according to inverse proportion of the numbers of cells for the remaining cell types.To achieve this, first integrated counts were exponentiated to make them non-log scale and slice_sample function from dplyr package (v1.0.7) was used with a replacement option.Then, these nonlog scale cell counts were summed up to create the pseudobulk profiles.For individual-based pseudobulk, AggregateExpression from Seurat package was used.6).AML tube 7 contains antibodies against CD41, CD25, CD34, CD117, CD42b, CD9, HLA-DR and CD45, but this tube has not been used to stain AML cases analyzed in this study.

Obtaining bulk transcriptomic data and their preprocessing
We have downloaded the non-normalized count matrices (htseq-counts) and the meta files of the four discovery cohorts (TCGA-LAML, BEAT-AML, TARGET-AML and TARGET-ALL) from https://portal.gdc.cancer.gov.For LEUCEGENE, count data was downloaded from their dedicated site (https://data.leucegene.iric.ca/)along with their provided meta data.All meta/count data were pre-processed using R (v4.1.0).For the meta data, genomic aberration labels were relabeled to the main AML WHO 2016 classes, non-AML samples were removed from the down-stream analyses, ELN-classes were relabeled according to ELN 2017 recommendations.For the count data, ERCC spike-ins and mitochondrial genes were removed, and the count matrix was then sorted according genes standard deviation in order to remove the duplicated genes that had less variation thus providing less information, and lastly the gene ensembl ids were converted to gene symbols.Before converting ensembl ids into gene symbols, the stemness score for each patient was calculated via count-per-million (cpm) normalized libraries, and these libraries were normalized using cpm function from edgeR package 33 (v.3.34.1).
Our 100 AML samples (LUMC) deposited to EGA with accession number EGAS00001003096 and they are accessible upon request.One hundred cryopreserved AML samples were selected from the Hematology Biobank of Leiden University Medical Center (LUMC) with approval by the institutional review board (no.B 18.047) and written informed consent were obtained according to the Declaration of Helsinki.QC benchmark analyses for these samples were done in our previous paper 11 .Therefore, we ran default HT-SEQ pipeline (v0.11.2) with paired-end option aligning fastq files to hg38 to obtain the count matrix.All above mentioned preprocessing steps (filtering, gene name conversion) were also conducted for these samples as well before deconvolution.

Deconvolution pipeline
To deconvolute the simulated pseudobulks and bulk RNA-seq AML patients, we used MuSiC 13 package (v0.1.1)as it benchmarked highly and consistently across different cell types at various settings 12 and had an easyto-use open-source (GPL-3) implementation (https://github.com/xuranw/MuSiC).We used non-log scaled count values as inputs and set the normalization option to false.Patients were assigned to the groups (e.g., GMP, CD14+ Monocyte etc.) according to their most abundant deconvoluted ECCs.In the heatmap (Fig. 2a), patients were re-ordered according to their ECCs within each assignment.

UMAP of estimated cell compositions
First, to obtain reproducible results with umap plots, we set a seed to 2 as UMAP procedure involves random initialization.Then, we ran umap function with default parameters from umap package (v0.2.7) and used the first two reduced dimensions to create the plots.All related figures were plotted using ggplot2 package (v3.3.5).

Drug resistance predictions via random forest
BEAT-AML has drug resistance data for 122 small-molecule inhibitors, we downloaded these from their manuscript (Supplementary Table S5).Then, each drug response was min-max normalized, then matched with their available RNA-seq samples.Next, drug resistances were predicted with the deconvoluted ECCs.Random forest algorithm from randomForest package (v.4.6-14) with default parameters was used for the predictions.Each sample within each drug was predicted at leave-one-out cross-validation settings.Then, for each drug, Spearman ρ values were calculated between predicted and actual drug resistance values.To stratify the drugs according to their primary diagnosis of WHO classification, samples within each diagnosis are selected and then each drugs predicted, and normalized drug resistance were associated, then the correlation and p-values were calculated using summary function in base R.

Venetoclax association analysis
First, each attribute was associated to standardized (min-max normalized) Venetoclax resistance from BEAT-AML study at univariate settings using lm function in R environment (v4.1.0)(Supplementary Table S12).Then, p-values and coefficients were calculated using summary function and then p-values were multiple hypothesis corrected using Benjamini-Hochberg procedure.For multivariate models, FAB classification, and ECC levels (except CD14+ Monocytes) were excluded as only 76 out of 460 samples of BEAT-AML had FAB classifications and as other ECCs are not independent of CD14+ Monocyte percentages (as the question is whether CD14+ Monocyte levels are independently predictive of Venetoclax resistance given BCL-2 expression in the same model).Again, p-values and coefficients were calculated with summary function and plotted in a volcano plot (Fig. 3f, Supplementary Table S13).
Survival analysis ggsurvival and survminer R packages was used for producing Kaplan-Meier curves.P-values were calculated with log-rank test.For ECC graph, patients were annotated with their most abundant cell type.Stemness score is split into low and high categories using the median value.

Proteomics sample preparation
Cell lysis, digestion and TMT labeling was performed as described in Paula et al. 34 .Cell lysis was performed using 5% SDS lysis buffer (100 mM Tris-HCl pH7.6) and 5 U benzonase nuclease (Thermo Scientific) with incubation at 95 °C for 4 min.Protein concentration was determined using Pierce BCA Gold protein assay (Thermo Fisher Scientific).100 µg protein of each sample was then reduced with 5 mM TCEP.Reduced disulfide bonds were alkylated using 15 mM iodoacetamide.Excess iodoacetamide was quenched using 10 mM DTT.Protein lysates were precipitated using chloroform/methanol; resulting pellets were re-solubilized in 40 mM HEPES pH 8.4 and digested using TPCK treated trypsin (1:12.5 enzyme/ protein ratio) overnight at 37 °C.Peptide concentration was then determined using Pierce BCA Gold protein assay.
The different samples, and reference samples, were arranged into five TMTpro 16plex sets.The peptides were labeled with TMTpro Label Reagents (Thermo Fisher Scientific) in a 1:4 ratio by mass (peptides/TMT reagents), total volume was 35 µL, for 1 h at RT. Excess TMT reagent was quenched with 5 µL 6% hydroxylamine for 15 min at RT.Samples corresponding to a TMT set were pooled and lyophilized.
Each TMT sample (80 ug) was fractionated by high pH reverse phase chromatrography on a Zorbax RRHD Eclipse Plus C18 2.1 × 150 mm 1.8micron column, at 800 ul/min using an Agilent1200 binary HPLC system, equipped with a UV detector.The mobile phases were 10 mM Ambic pH 8.4 (A) and 10 mM Ambic/Acetonitrile 20/80 pH 8.4 (B).The gradient was from 2% to 90%B in 35 min.20 fractions were collected in a circular fashion, i.e., collection per vial for 20 sec before moving to the next collection vial.After collection in the last vial collection is continued in the first vial.Fractions were subsequently freeze dried.

Mass spectrometry
TMT-labeled peptide fractions were dissolved in water/formic acid (100/ 0.1 v/v) and analyzed by on-line C18 nanoHPLC MS/MS with a system consisting of an Ultimate3000nano gradient HPLC system (Thermo, Bremen, Germany), and an Exploris480 mass spectrometer (Thermo) as in Rossi et al. 35 .Fractions were injected onto a cartridge precolumn (300 μm × 5 mm, C18 PepMap, 5 um, 100 A), and eluted via a homemade analytical nano-HPLC column (50 cm × 75 μm; Reprosil-Pur C18-AQ 1.9 um, 120 A (Dr. Maisch, Ammerbuch, Germany)).Solvent A was water/ formic acid 100/0.1 (v/v).The gradient was run from 2% to 40% solvent B (20/80/0.1 water/acetonitrile/formic acid (FA) v/v) in 120 min.The nano-HPLC column was drawn to a tip of ∼ 10 μm and acted as the electrospray needle of the MS source.The mass spectrometer was operated in datadependent MS/MS mode with a cycle time of 3 s, with the HCD collision energy at 36 V and recording of the MS2 spectrum in the orbitrap, with a quadrupole isolation width of 1.2 Da.In the master scan (MS1) the resolution was 120,000, the scan range 350-1200, at standard AGC target @maximum fill time of 50 ms.A lock mass correction on the background ion m/z = 445.12003was used.Precursors were dynamically excluded after n = 1 with an exclusion duration of 45 s, and with a precursor range of 20 ppm.Charge states 2-5 were included.For MS2 the first mass was set to 110 Da, and the MS2 scan resolution was 45,000 at an AGC target of 200% @maximum fill time of 60 ms.

Proteomics data processing and down-stream analysis
In a post-analysis process, raw data were first converted to peak lists using Proteome Discoverer version 2.4 (Thermo Electron), and submitted to the Uniprot database (Homo sapiens, 20596 entries), using Mascot v. 2.2.07 (www.matrixscience.com)for protein identification.Mascot searches were performed with 10 ppm and 0.02 Da deviation for precursor and fragment mass, respectively, and the enzyme trypsin was specified.Up to two missed cleavages were allowed.Methionine oxidation and acetyl on protein N-terminus were set as variable modifications.Carbamidomethyl on Cys and TMTpro on Lys and N-terminus were set as fixed modifications.Protein and peptide FDR were set to 1%.Normalization was on the total peptide amount.The 5 TMT-16plex analyses were normalized to each other by the bridge samples.
First, the abundance data is log-cpm transformed to stabilize variance among samples and then to ensure dealing with the technical batch effects, we ran removeBatchEffect function from limma package (v3.48.3) providing batch information.Then, we overlaid the transformed protein abundances onto a PCA plot to ensure that there is no batch related clustering, as can be observed from the positioning of the TMT control samples before (Fig. S6A) and after correction (Fig. S6B).The down-stream analyses were then done using only the primary AML (n = 39).Next, the transformed BCL-2 protein abundance was associated with the BCL-2 expression from LUMC cohort (Fig. 4a) and CD14+ Monocyte percentage within each sample (Fig. 4b) (Supplementary Table S14).R 2 and p-values were calculated using stat_-poly_eq function from ggpmisc (v0.4.5) package and lines were drawn with geom_smooth function from ggplot2 package (v3.3.5)via method option set to linear model.

d
Deconvolution of in-silico pseudobulks from Human Cell Atlas (HCA) CD14+ Monocyte percentage with flow cytometry and ECC Comparison of flow cytometry and estimated cell composition ECCs of TARGET-AML and - Heatmap of deconvolution results for all patients UMAP plot of cell type assignments FAB classifications agree with the estimated cell compositions (ECCs)

Fig. 2 |R 2 R 2 Fig. 3 |
Fig. 2 | Capturing differentiation stages with ECCs. a Heatmap showing the deconvoluted fractions for of five AML cohorts (TCGA-LAML, BEAT-AML, TARGET-AML, LEUCEGENE and LUMC).Each column represents a sample, which is deconvoluted into 22 cell types (bottom part with blue) via healthy single cell BM reference.Top part of the annotations shows the provided meta data (FAB, WHO and ELN classifications, origin, sex and reported blast percentages) for these cohorts.Each patient was assigned to its most abundant deconvoluted cell type (AML phenotype) and samples within each assignment were sorted according to the assigned phenotype.b UMAP plot of deconvolution levels annotated via most abundant cell type as the left panel; and on the right, ECCs of 6 out of 22 cell types were shown in continuous scale (HSC, CD14+ Monocytes, Late Erythrocytes, GMP, EMP and cDC).c Similar to UMAP plot in (b) but colored with FAB classifications (left panel), and the fractions of cell types for each of these FAB classifications (right panel).
10L cases were stained with fluorescent antibodies and analyzed by flow cytometry for diagnosis, prognosis, and disease monitoring of AML in the diagnostic laboratory of the department of Hematology in the Leiden University Medical Center.The flow cytometric test has been developed and performed according to EuroFlow standard operating procedures (www.euroflow.org)10.EuroFlow antibody combinations have been tested against references databases of normal cells from healthy individuals and allow multidimensional identification and distinction of aberrant cells from normal cell populations.The flow cytometric test includes 8 tubes with different antibody combinations, i.e. one ALOT tube (acute leukemia orientation tube) and 7 AML tubes (AML1-7).